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1 Introduction 

Over the past three years, the Computational Mechanics Co., Inc. has designed and produced 
a new computational tool for modeling internal flows in solid rocket motors as part of the 
project “Internal Flow Analysis by Highly-Accurate Adaptive Computational Methods” 
(NAS8-37682). During the course of this effort several new algorithms for studying the 
effects of moving boundaries on flow characteristics within solid propellant rocket motors have 
been developed and tested on both two-dimensional planar and axisymmetnc computational 
domains. Of particular interest has been the development of a receding boundary algorithm 
which successfully models the changing flow domain during the erosion of the propellant 
within a solid rocket motor. This method relies on adaptive finite element methods that 
combine node relocation techniques (r-methods), mesh refinement techniques (h- methods), 
and moving boundaries to form a very powerful analysis package. 

In addition to the research and development efforts in the area of moving boundaries, 
considerable effort has also been focused on the following topics: 

• Implementation of a fully adaptive implicit/explicit finite element methodology to op- 
timize the computational effort required to advance the solution forward in time. 

• Formulation and implementation of a generalized algebraic turbulence model and data 
structure compatible with adaptive methodologies and applicable to completely un- 
structured grids. 

• Testing and validation of the computational algorithms for several benchmark prob- 
lems. 

The success of the two-dimensional and axisymmetric analysis package in modeling the 
benchmark problems suggests that extensions to realistic three-dimensional flows in cavities 
with eroding boundaries is both feasible and a natural extension of the project currently 
underway. The final section of this report describes some possible extensions of this project 
in terms of enhancements of the two-dimensional code and the development of a fully three- 
dimensional code for use in the design and analysis of solid rocket motors. 

The remainder of this report presents a detailed description of the theoretical formulation, 
numerical methods, and representative examples of the analysis of complex flow phenomena 
occurring in solid rocket motors. The equations that model flow phenomena for this class 
of problems are derived in Section 2. The associated boundary conditions are then briefly 
discussed in Section 3. An implicit/explicit flow algorithm is presented next in Section 4, 
while the moving boundary and remeshing algorithms are highlighted in Section 5. Section 
6 introduces adaptive strategies and the related data structures. The implementation of a 
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simple algebraic turbulence model is presented in Section 7 along with a discussion of the 
data structure and storage requirements. Following the discussion of the algebraic turbulence 
model are several numerical examples which demonstrate some of the modeling capabilities 
available in the code. 


2 Governing Equations 

For the class of the fluid dynamic problems we are going to solve, the computational domain 
may be continuously changed due to boundary motion. It is well known that for fluid 
dynamic problems involving grid motion, the original Navier- Stokes equations derived from 
the Eulerian approach need to be modified by using the Arbitrary Lagrangian-Eulerian 
(ALE) formulation, see reference [ 1 ] for detailed derivation. It should be noted that the 
ALE formulation derived in [ 1 ] is expressed in the form of Cartesian tensor notations. To 
derive its axisymmetric counterpart, we first convert the Cartesian tensor notations to vector 
operator forms. These vector operators can be readily transformed into axisymmetric form 
using orthogonal curvilinear transformations. 


Two-Dimensional Formulation of the Navier-Stokes Equations 


The time-dependent Navier-Stokes equations, including grid motion, can be expressed in 

tensor notations as _ 

Bp , Sri*—,) =0 


dt + 

dpui d 
dt + dxj 
de d 
dt ^ dx{ 


d x, 

pui(u 3 - uf) - r.j + opui = 0 
e(ui - uf) - UjTij + q,] + ere = 0 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


where p denotes the density of the fluids, u, is the i-th component of flow velocity, uf is the 
i-th component of grid velocity, and e represents the total energy per unit mass. The shear 
stress tensor r,j can be expressed as 

. , r duk ( du{ duj\ 
r tJ = -fa + XS l3 — + p 


and the heat flux < 7 , is 


9. = 


-k 


dT 

dxi 


Here p and A represent the molecular and secondary viscosity, respectively. The source terms 
in equations (2.1 )— (2.3) are due to the grid volumetric dilatation, a = V u G . For simplicity 
and clarity, we will omit these terms in subsequent discussions. 
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Axisymmetric Formulation of the Navier-Stokes Equations 

In vector operator form, equations (2. 1-2.3) can be expressed as 

fy + v . p( v- V°) = 0 


(2.4) 


+ V • ( V — V G ) ® pV + Vp — V(2fi + A)V ■ V 
3‘ (2.5) 

+ V x /x(V x V) = 0 

— + V • e(V - V G ) + V pV - V • (AW ■ V) 
dt ( 2 . 6 ) 

- V ■ fiV(V -V) + V - {nV xV xV) + V - q = 0 

Here V is the velocity vector and V G represents grid motion. 

In orthogonal curvilinear coordinate systems (ei, 62,63), these vector operators can be 


■^(h 2 h 3 v i) + JL(h x h 3 v 2 ) + ^-{h,h 2 v 3 ) 
axi ox 2 0x3 J 

h\e\ h 2 e 2 h 3 e 3 

d d d 

dxi dx 2 dx 3 

hiV-i h 2 V 2 h 3 V 3 

where /i, are metric scale factors. 

For a cylindrical coordinate system, 

h x = 1, h 2 = r, h 3 - 1 

By expanding these vector operators in equations (2.4-2.6), and neglecting all terms 
related to 6, the axisymmetric form of Navier-Stokes equations is obtained: 


expressed as: 

V 


e, 


V • V = 


\_d_ 

hi dx x 

1 

hi h 2 h 3 


V x V = 


1 

h\h 2 h 3 
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Equation of Continuity 


dp_ d[pUr) djpUz) /ntr _ Q 

dt dr dz r 


Multiplying equation (2.7) by r and rearranging terms, we have 


djpr) djpru^) d{pru 2 ) 

dt dr dz 


Equation of motion (r-component) 



where 


r rr = 2 M ^ + AV-V, 
( du r du z \ 

T " = " Ur + ir 


dr dz r 


X = — 2/r/3, if Stokes’ hypothesis is applied 






Multiplying equation (2.9) by r and rearranging terms, we obtain 


where 


and 



( 2 . 10 ) 


Equation of motion (z-component) 



where 


Tzz 


2^ + AV • V 

az 


Tzr 



du z 

dr 


) 
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Multiplying by r and rearranging terms, we have 



( 2 . 12 ) 


Energy Equation 



Multiplying by r and rearranging terms, we have 



To make these equations consistent with the two-dimensional planar counterpart, the 
notation is changed from 2 to x and r to y. Depending upon the arrangement of the terms 
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in the governing equations, two different forms may be obtained, 
from equations (2.7), (2.9), (2.11), and (2.13) as 


ii + FI - FI + (S c -S v ) = 0 


where 


u = [p pui pu 2 e] T 
p (u, - up) 

pui ( Ui - Up] + p6u 

pu 2 ( Ui - u G ) + pS 2 i 

£ (ui ~ Up) + UiP 


F v , = 


0 

T 1, 



P (u2 ~ ^P) 
pu^ (u 2 - up) 
pu 2 (u 2 - up) 

£ (u 2 — Up) + pU 2 



and the viscous shear stresses are 


0 

T 12 

2, - fj 


Tn = + Au 2|2 


The first form is concluded 


(2.15) 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 
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T 22 — 

PRU 2,2 + Au 1; i 

(2.22) 

r 21 = 

T 12 = P (Ui, 2 + U 2 ,l) 

(2.23) 

PR = 

(2 p + A) = longitudinal viscosity 

(2.24) 


It is observed that in this version, both convective and viscous fluxes are identical to 
the two-dimensional planar counterpart and the additional source terms include the grid 
velocity. The advantage of applying these equations in the code development is that only 
the source terms need to be incorporated into the code and, therefore, the code integrity 
and efficiency may be maintained. Unfortunately, it can be shown that this form is not 
conservative, and thus unacceptable for the numerical implementation. 

The second form is derived from equations (2.8), (2.10), (2.12), and (2.14) as 


Z + Fl-Fli + S'-S^ 0 


(2.25) 


where 


2 

u = yu = y[p pu-i pu 2 e] 


F, = yF\ = y 


p (u,- -up) 

pu-i (uj — Up) + p6u 

pu 2 (u; — Up ) + pS 2 i 
[ £ (ui - up) + u { p 


F x = yF) = y 


0 

Tli 

T* 


S c = — [0 0 p 0] 

s v = - 


o 0 + A u m , m 0 

y J 


(2.26) 


(2.27) 


(2.28) 


(2.29) 

(2.30) 
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and the viscous shear stresses are 


Tn 

A u 2 

= Pru 1,1 + Au 2 ,2 H 

y 

(2.31) 

r 2 2 

A u 2 

= PRU 2,2 + Alii 1 H 

y 

(2.32) 

T21 

= Ti2 = P (ui,2 + «2,l) 

(2.33) 

Pr 

= (2 p. + A) = longitudinal viscosity 

(2.34) 


This form has the advantage that the grid velocity makes no contribution to the source 
terms and there is no formal difference between the axisymmetric and planar formulations 
due to the grid velocity. Although this form is conservative and will be employed for the 
code development, the finite element interpolation based on u = y (p pu\ pii? £) will cause 
severe numerical problems near the axis of singularity, y = 0, since the conservation variables 
u-(p pui pu 2 e) T for all nodes that lie on this axis can not be recovered directly from u. 
We have resolved this problem by interpolating u and y separately, see Section 4.1 for a 
detailed description. 

3 Boundary Conditions 

It is well known that the well-posedness of fluid dynamic problems can not be established 
by solely considering the governing equations without investigating the associated boundary 
conditions [2, 3, 4, 5]. Moreover, the accuracy and the rate of convergence of a numerical 
algorithm are also affected by the proper treatment of the numerical boundary conditions 
[6-10]. On the other hand, the ability of a CFD code to simulate fluid dynamic problems 
is strongly dependent on the flexibility of treating various types of boundary conditions. In 
this report, instead of presenting a detailed mathematical background about these boundary 
conditions, we will rather point out various types of boundary conditions that are most 
often encountered in fluid dynamic problems and discuss how these boundary conditions are 
implemented within the context of finite element methods. 

The following boundary conditions are common to most fluid dynamic problems: 

1. Open boundaries or inflow/outflow boundaries. 

2. No-flow or no-penetration boundary. 

3. No-slip isothermal/ adiabatic boundary. 

4. Porous wall boundary with mass injection. 
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5. Moving boundary with prescribed velocity. 

Boundary condition type 1 is usually encountered when artificial boundaries are intro- 
duced to a problem in order to reduce the computational domain. In this case, boundary 
conditions have to be treated carefully to handle wave reflection and disturbance propagation 
for both accuracy and rate of convergence considerations [6-10]. 

Boundary condition type 2 is popular for inviscid flow problems and for problems involv- 
ing symmetric geometry. On this type of boundary, flow is prohibited from penetrating the 
boundary and is allowed only in the direction tangent to the boundary. This condition is 
quite natural for solid boundaries immersed in an inviscid flow. 

Boundary condition type 3 is a typical boundary condition for viscous flow problems. 
The isothermal boundary simply implies that the temperature of the wall is a constant, 
while adiabatic boundary means that the heat flux across the wall is zero. 

Boundary conditions 4 and 5 are special features associated with the numerical simulation 
of eroding or moving boundaries. Generally speaking, the eroding surface of a solid propellant 
may be regarded as a porous wall with boundary motion. The burning rate and the mass 
flow rate across the surface are, theoretically, dependent on the local flow conditions such as 
pressure and temperature and the composition of the solid propellant. A notable empirical 
formula known as Saint-Robert’s law which correlates the burning rate and chamber pressure 
is given by [11]. In cases where chemical reactions and combustion phenomena are not 
considered, the surface of the solid propellant is usually treated as a porous wall with mass 
injection, the mass injection rate being obtained from experimental data. 

To avoid ambiguity in identifying the type of boundary condition at a node where two 
different types of boundaries intersect, a hierarchical procedure is set up as follows based on 
their relative priority: 

1. porous wall boundary with boundary motion, 

2. no-slip boundary condition, 

3. no-flow or no-penetration boundary condition, and 

4. open boundary condition. 

A detailed discussion of the theoretical formulation and implementation of most of these 
boundary conditions can be found in reference [20]. In this work we will discuss boundary 
conditions specific to the solid rocket booster applications, namely porous/burning wall 
conditions (Section 4.4). 
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4 Numerical Methods 

Over the past few years, significant progress has been made in the development of new 
computational methods for solving compressible Navier- Stokes equations. The approaches 
described in the literature vary from fully explicit algorithms, which are computationally 
inexpensive but often severely limited by stability restrictions, to fully implicit algorithms, 
which are unconditionally stable but are much more expensive per time step. The selection 
of which type of algorithms is optimal for a given application is generally not known a prion 
and may in fact change as the features of the flowfield develop. As a result, a domain 
decomposition approach is gaining popularity in the CFD community which employs an 
explicit formulation in one region of the mesh and a fully implicit formulation in another, 
see references [12-15,20]. 

In this section, a family of implicit/explicit finite element algorithms developed by 
Tworzydlo, Oden, and Thornton [20] will be generalized for the solution of axisymmetric 
problems with moving meshes and porous/burning wall boundary conditions. Within this 
family, a fully explicit method or various versions of implicit algorithms can be obtained by 
appropriate selection of implicitness parameters. This general finite element algorithm is 
combined with adaptive mesh refinement (as presented in references [18,19]) and adaptive 
selection of implicit/explicit zones within the computational domain. Several approaches to 
the selection of implicit and explicit zones are presented. 

Following the description of implicit/explicit schemes, the artificial dissipation added to 
the system to suppress possible spurious solutions (due to shocks) will be outlined. The final 
solution describes the implementation of the porous wall boundary condition and pressure 
outflow boundary condition. 


4.1 A General Family of Implicit Taylor-Galerkin Methods 

In this section a general family of implicit Taylor-Galerkin methods will be derived. This 
family is based on a combination of second— order Taylor series expansions in time with a 
Galerkin approximation in space (one-, two-, or three-dimensional). Several implicitness 
parameters are introduced, so that, depending on the particular choice, a fully explicit scheme 
or a variety of implicit schemes can be recovered. The family of algorithms presented here 
is a generalization to axisymmetric problems of the algorithms developed previously. 

Second-Order Taylor Expansion in Time 

Assume that the solution u n is given at the time moment t n and the solution at time t n+l is to 
be calculated. Formally, the values of the solution at moments t n and f n+1 can be expressed 
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by the second— order Taylor expansion around an arbitrary moment t + (see Fig. 4.1), where 
a is the implicitness parameter with values between zero and one: 


u n+1 = u n+ ° + (1 - q)A tu nJra + (1 - a) 2 ^-w n+a + 0{At 3 ) 


U n = u n+a - aAti n+ ° + a 2 ^-u +a + 0(Af 3 ) 


(4.1) 


By subtracting these two formulas one obtains a formula for the increment of the solution 
between steps n and n + 1: 


Afi = u n+1 - u n — Aw n+ ° + (1 — 2 a) 


A f 2 ^n+Q 

-u + O(Ar) 


(4.2) 


Now it is easy to observe that: 


n+a = ^ + _ /0 ) Af ) 


so that — still preserving second-order accuracy — a second implicitness parameter 8 can be 
introduced into equation (4.2): 


Am = Atu" + ° + (1 — 2a)— u 0 + 0(At 3 ) (4-3) 

The next step of the derivation is to express the quantities evaluated at time moments t 
and t n+ @ by quantities evaluated at the basic steps t n and f n+1 . It is easy to show, using a 
Taylor series expansion, that: 

u n+Q = u" + aAu + 0(Af 2 ) 

u +p = u +8Aii + 0(At 2 ) d 

Substituting these formulas into equation (4.2) yields a two-parameter expansion: 


Au = At (tr" + aAu) + (1 — 2a)—— (u + fiAu) -f 0(Af 3 ) 


(4.4) 
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Now, following the original idea of Lax and Wendroff, the original equation (4.1) will be sub- 
stituted to equation (4.4) to replace time derivatives by space derivatives. This substitution 
yields a formula for the first derivatives: 


£ = Fl, - Fl + S" - S c 


(4.5) 


AS = AF*^ - AF- t - + AS V - AS C 

= (%, AS.. j) . + (P.Aii) . - (AA«) . 
+ QjAu,j + TAu - BAS 
and for the second derivatives: 

j-v ■ v ■ c 

u = F iti - F,, + S-S 

= (RiAj + Piu) . - , 

+ (QjUj + Tu) + Bu 

where the relevant Jacobians are defined as: 

Rij - 

A t = 


dF- 

D 

dF > 

diij 

1 * i 

du 

dF- 

du 

> Q i 

dS v 

duj 

dS v 


dS c 

du 

> & 

du 


(4.6) 


(4.7) 


T = 


It can be shown that the multiplier, y , will not affect the linearization of nonlinear invis- 
cid fluxes. Moreover, Jacobians J3,j will remain the same as its two-dimensional planar 
counterpart. This implies that 

Bjj = B,j ) A-j A t 

where Rij and A, are the planar counterpart and can be found in reference [23,24]. The 
Jacobians Pi can be expressed as the sum of two components as 


P, = Pi + Pf 
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where Pi is the two-dimensional planar counterpart and the additional term -P, is due to 
the axisymmetric formulation. Details of P;, Q : , T, and B are given in Appendix A. 

After further substitutions, the second order time derivative can be expressed as 

fi = [a, (F” * - r u + s- - S') .] . 

- -f m + s--s‘)] . 


+ 


(4.8) 


+ [« ; (?M-Fb + 

- [r(F* w -Fl. t + S v -S')l 

- B (F kik - T iit + S* - 5')] 

Since bilinear elements are used in the finite element code, any spatial derivatives higher 
than order 3 will be neglected in equation (4.8). This results in the following approximation 
for the second order time derivatives. 

£ = (aKJ, + (*$%• + BF c kk + BS c + 0(n,k) 

= ( AiA k u ik ) x + (A,Bw). + BA k u <k + BBu + 0{n,k) 

and the incremental form of the second order time derivative can be expressed as 
Af2 = + (A,BAu) . + BA k Au tk 

+ BBAu 4* 0(Ai) + o(/i, k ) 

Substituting (4.5)-(4.6) and (4.9)-(4.10) into (4.2) and regrouping explicit and implicit 
terms, we get 


(4.9) 


(4.10) 


Au + aAt ( AiAu ) • + aAtBAu - 7A t ( RijAuj ) t 
- 7 At (P, Au) - ~fAtQ : Auj - lAtTAu 


(1 -2a)0At 


- (A t A k AU' k ) t + ( A t BAu ) . 


+ BA k Av,' k + BBAu] = At (F ^ - F iti + S v - S c ) 

+ (1 ~ 2 , a ^ [(Ai?;.,) + (A, S'), + + BS'] 


(4.11) 


+ (1 - 2 a)0(n, k)At 2 + (a - 7 )0(/i, k)At 2 + 0(At 3 ) 
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Notice that a third implicit parameter, 7, was introduced in (4.11) to control the implicitness 
of the viscous terms without affecting the order of accuracy in time. 

Variational Formulation 

Neglecting higher order terms and multiplying by an arbitrary test function v through (4.11) 
and integrating over the domain we get 

/ (Ait • v + aAf[(-A t Au)v, t + (BAu)v) 

J n 

+ 'yAt[{R ij Auj)v ii + ( P,Au)v , - (QjAuj)v 
- (TAu)m] + (1 - 2a)0A L [( i 4. i 4 ifc Att,fc)i>, t - 

+ (A t B Au)t7,i - (BAj Au tj )v - (BBAu)v}}dtt 

+ J^{aAt(riiAiAu)v - 'yAt[(n x R l] Au iJ )v (4.12) 

+ (n,P,Aw)f] - — ~ 2 ^ A< {( n ,A t BAu)v + (n 1 A 1 A J Au J )M]}ds 

= / {At[(F- - F V i)v,i - (S c - S u )v] 

Jn 

+ — — ^ — [-{AiAjU t j)v t i - (A,s>, t + (BAjUj)v 

+ (BS>]}<fft + |Atn,(F- - ?■) v + (1 ~ 2 -^^- {{n l A i A j u <3 )v + (n, A,S>] j ds 

In the variational formulation (4.12), the unknowns selected for finite element interpolation 
is u = yu, or in discrete form 

j v 

u(x) = ^2/^/(x) (4.13) 

/=i 

The choice of u for the finite interpolation leads to the following computational problems: 
u becomes zero as y —> 0. This implies that the conservation variable u becomes singular, 

u = lim — = lim — (4.14) 

y -> o y y—*0 y 
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To avoid this singularity problem, we will interpolate y and u separately for u so that 
u becomes the actual vector of unknowns in the problem. The resulting variational form 
becomes 


f {yAu ■ v + aAf[y(-A"Au)r,- + y(BAu)v) 

In 

+ 'yAt{y{R ij Au tJ )v ti + {Ri : yj Au)v^ + y(P t Au)v it 
- yiQjAu^v - (Qjy t jAu)v - y{TAu)v] 

+ — — ^Y^- iyiAiAjAu^v^ + (AiAjyjAufa'i 
+ y(AiBAu)v, t - y(BAjAuj)v - (BA } y,jAu)v - y(BBAu)v]}dCl 


f {aAty(rnA t Au)v - 7At[y(n,H lj Au J )r + (n.f^t/j Au)r + y(n,P t Au)v ] (4.15) 

JdQ 

— — — ^y^ n .^.A^Auj)v + (riiAi Ajy tJ Aw) • v + y(n t AiBAu)v]ds 
= J ' A t[y(F C t - F“)v,i - ( S c - S v )v] - (1 - 2 2 ° )At [y{A,A 3 u 3 )v + (. A t Ajy yj u)v ti 


+ AiS c v ti - y(BAjUj)v - ( BAjy tJ u)v - ( BS c )v]dSl 

+ / {Aiyn,-(F* - F‘)v + — — 2 r ^ At [y(niAiAjUj)v 

JdQ J 

+ UiAiAjyju) ■ v + (mAiS c ) ■ r]}ds 

Since matrices B, T 1 , Q, and P contain multiplier j and/or ^ (see Appendix A) and i / ,1 — 0, 
y 2 = 1, the integrals containing these terms in equation (4.15) can be regrouped as: 
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1. Left hand side interior integrals 


- J f [6 (S 2 + BA 2 ) +c(Q 2 + T 2 )\ 

+ J (aB — cTi)Au ■ vd£l + J I(yAu • v)dQ, 

+ / (-aAi + cPi)(yAu-V'i)d£l 

Jn 

+ J [(cR,^ + bAiA 2 ) + (bA x B) + cP ?} (Au ■ v ,;)<tf) 

- J (cQj + bBA^j {Auj • v) dfi, 

+ l (cRij + bAiAj) (yAujV'i) dQ 
Jn 

A , (1 — 2a)/3At 2 . 4 

where a = aAt , b = — , c — jAt 


2. Left hand side boundary integrals 


/ (an t A{ - cn t P t ) (y Au ■ v)ds 
Jd n 

- JJ cn t Ri 2 + bniAiA 2 ) + bniAiB - cn.P^] (Aw • v)ds 



+ bn,AiAj) (yAujv) ds 


(4.16) 


(4.17) 
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3. Right hand side interior integrals 

f At {S\ - S c ) ■ vdSl 

Jn 

+ J f |aiS> + — — 2 °^ A< [( BA, ) u ■ V + BS C • »]| 

+ J a & 

- J a | At (F”' 1 • «,) + ^ ~ 2 °) AI (a, (4 2 u + S') ■ i>,| | a 

+ / ( ‘ ~ 2a) — BXjK-.IJI 
Jo 2 

- I ^—^Y^- AiAj (ttj • »,,•) 


where 5j = — [0 0 Au^* 0] r 

S 2 = -[0 0 /ifiU 2 0] T 

F*" 1 is the Cartesian counterpart of viscous fluxes and F* A is the additional 
fluxes associated with the multiplier, A, where 

F\ A = [0 Xv 0 Auu] T 

F v 2 a = [0 0 Au Au 2 ] T 


4. Right hand side boundary integrals 

[ At [ym (*7 C - F c t ) • v] ds + / A tn x F vA ■ vds 
Jdti 

+ [ — ~ (niAiA 2 ) (u • t>) ds 

JdQ 2 

+ f AzAB^ niAi (s‘-v)is 

JdQ 2 

+ / — — ‘-^—— n x A x A J (; yuj • ti)ds 

JdQ 2 


(4.18) 


(4.19) 
viscous 

(4.20) 


(4.21) 
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4.2 Implicit/Explicit Procedures 

The basic idea of implicit/explicit algorithms is simple: combine the two methods to take 
advantage of the superior features of each. The major advantage of the explicit method 
is that element computations are relatively inexpensive and simple. Unfortunately this 
method suffers from stability limitations of the time step, which in some problems leads to 
prohibitively large numbers of time steps. 

The implicit algorithm allows for an application of larger time steps than the explicit 
method. Moreover, due to the existence of implicit boundary terms, it offers easy and 
straightforward control of natural boundary conditions, particularly those involving the vis- 
cous fluxes. An additional advantage is that with larger time steps no explicit artificial 
dissipation is necessary, which is very important in the calculation of boundary fluxes, par- 
ticularly by wall heating rates. The major disadvantage of implicit methods is a much higher 
cost of element operations and a more complex and expensive solution of the resulting system 
of equations. 

In this section, the formulation and numerical implementation of an adaptive implicit/ 
explicit algorithm for compressible flows is presented. The algorithm will be based on the 
general family of Taylor- Galerkin methods discussed in the previous sections. 


Formulation of Implicit/Explicit Schemes 

The algorithms for determining the partition must be designed so as to preserve stability, the 
conservative properties, and the required order of approximation. We begin by partitioning 
the domain f 2 into subdomains Cl^ and f where explicit and implicit schemes are to be 
applied, respectively. 

Some of the possible procedures are examined below. 

Procedure I 

The first possible approach, applied for example by Hassan, Morgan, and Peraire [25], is 
based on the following two steps: 

1. Perform the explicit step computations on all nodes in the mesh (0^ = fl). 

2. Perform the implicit computations in the subregions, where the stability criterion for 
the explicit scheme: C < 1 is violated. The solution in remaining nodes is “frozen” at 
this step. 
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This simple procedure has one basic disadvantage: it appears to be nonconservative and 
this disturbs the regularity of the solution in the transition zone. This is caused by the 
fact that during Step 2 the “frozen” explicit nodes impose the actual Dirichlet boundary 
condition on the edge of the implicit zone. Prescribing the Dirichlet boundary condition 
for the Navier-Stokes equation means that there must exist an (external) source of fluxes 
to support the prescribed state of the solution. Since no such external source exists within 
the computational domain, the solution will not be conservative across the implicit/exphcit 
line. Due to enforcement of the Dirichlet boundary conditions, the solution may exhibit a 
“ramp” or “kink” along this line. 


Procedure II 


Again fi is considered to be the union of two subregions and 0^ ^ (see Fig. 4.2), such 
that: 

n {E) n q (/) = t E i , n {E) U 0 (/) = D 

It is convenient to assume that the interface between the two regions coincides with the 
element boundaries. 

It can easily be observed that the differential equations to be solved on the two subregions 
are different due to different implicitness parameters applied in each zone: ^ 

the implicit zone and = 0 in the explicit zone. Therefore, the variational 

formulation (4.15), based on the assumption of constant implicitness parameters, cannot 
be applied to the domain 0. Instead, it can be applied separately to each subdomain with 
additional continuation conditions across the interface. These conditions represent continuity 
of the solution and satisfaction of the conservation laws across the interface and are of the 
form: 



= 

p(E)C 

= f\ I)C 

a (E) 

= A™ 

F (E)V 

= Fl! )v 


on r^/ 


(4.22) 


where index n refers to the outward normal for the corresonding region (n (E) = -n (/) ). The 
continuity requirement also pertains to the test function, so that — v. Note 

that for general weak solutions of Euler equations the solution v. need not be continuous 
across the interface. However, for regularized problems and finite element interpolation, the 
continuity of u is actually satisfied. 

If the variational statement is formulated for this problem, then in addition to interior 
integrals for each subdomain and regular boundary integrals, jump integrals across the in- 
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terface appear in the formulation. Note that in the practical implementation of this scheme 
we set all the interface components to zero. This procedure preserves the continuity of fluxes 
and time accuracy across the interface up to the first order. 


Procedure III 

The last procedure considered here is based on a generalization of the weak formulation, 
according to which the implicitness parameters are not constant, but are continuous functions 
of the position x. For the finite element computations, it is convenient to limit the choice of 
these parameters to the finite element subspace, so that: 

N 

a{x) = a /'M x ) 

1=1 

and the same holds for the other implicitness parameters 0,7 and 8. With this assumption 
additional terms show up in the variational formulation. Note that there are no additional 
boundary terms resulting from the variable implicitness. 

The above approach seems to be the most general and clear, with no ambiguities con- 
cerning the interface conditions. To make the practical application cheaper, the implicitness 
coefficients are held constant in most of the elements, and the additional terms are actually 
evaluated only in the transition zones. 

Selection of Implicit and Explicit Zones 

The basic criterion for selection of implicit and explicit zones is simple: for a given time step 
all nodes which violate the stability criterion for an explicit scheme should be treated with 
the implicit scheme. According to this criterion, several options for an automatic adaptive 
selection of implicit/explicit zones were implemented: 

1. User-prescribed time step DT : 

Within this option, the user prescribes the time step. All nodes satisfying stability 
criterion for the explicit scheme (with s certain safety factor) are explicit. This means 
that all the elements connected to these nodes are treated with the explicit scheme. 
On all other elements the implicit scheme is applied. 

2. Prescribed maximum CFL number: 

In this option, the user prescribes the maximum CFL number that can occur for 
elements in Q. The time step is automatically selected as the maximum step satisfying 
this condition. The choice of a maximum CFL number may be suggested by the time 
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accuracy arguments or the quality of results (it is known that, for a Taylor-Galerkin 
scheme, too large a CFL number tends to smear shocks). 

3. A prescribed percentage of the domain is implicit. 

In this version, the user specifies the fraction of the domain which is to be treated im- 
plicitly. The elements with the strongest stability limitation (usually the smallest ones) 
are treated implicitly, the others are explicit. The time step is selected to guarantee 
stability of the explicit zone. 

4. Minimization of the cost of computations. 

In this option, the time step and the implicit /explicit subzones are selected to minimize 
the cost of advancing the solution in time (say one time unit). The algorithm is based 
on the fact that, for an increased time step, an increasing number of elements must be 
analyzed with the (expensive) implicit algorithm. The typical situation is presented in 
Fig. 4.3, which shows for different time steps the relative number of nodes that must 
be treated with the implicit scheme (to preserve stability). On the abscissa, the A t F £ 
denotes the longest time step allowable for the fully explicit scheme (with certain safety 
factors). A t FI denotes the shortest time step requiring a fully implicit procedure. The 
relative number of implicit nodes increases as a step function from zero for At < A t F E 
to one for At > t F j. Now assume that the ratio r of the computational cost of 
processing one implicit node to the cost of processing one explicit node is given. This 
ratio can be estimated relatively well by comparing the calculation time of element 
matrices and adding, for implicit nodes, a correction for the solution of the system 
of equations. Then the reduction of the cost of advancing the solution in time with 
the implicit/explicit scheme, as compared to the fully explicit scheme, is given by the 
formula: 

R{At) = ^£^ n (E) +rn (^ 

Typical plots of the function R{At) are presented in Fig. 4.4. Shown here are the two 
cases: 

(a) the case of a small difference between fully explicit and fully implicit time steps 
an almost uniform mesh 

(b) the case of a large difference between fully explicit and fully implicit time steps 

Note that in either case, restrictions on the length of the time step should be applied, 
for example, from the maximum CFL condition. Otherwise the cheapest procedure 
would always be to reach the final time with one implicit step. 

From the plots in Fig. 4.4, the following observation can be made: for an essentially 
uniform mesh, the mixed implicit/explicit procedure does not provide savings of the 
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AT e - fully explicit time step 



Figure 4.3: Relative number of implicit elements for increasing time step. 
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computational cost — either a fully implicit or fully explicit scheme is the cheapest de- 
pending on the time step restriction. On the other hand, for very diverse mesh sizes the 
mixed procedure provides considerable savings. This means that the effectiveness of 
the mixed implicit/explicit scheme will be the best for large-scale computations with 
both very large and very small elements present in the domain. In the practical imple- 
mentation of this method, the approximation of the function R{At) is automatically 
estimated for a given mesh. Then, the time step corresponding to the smallest R(At) 
is selected automatically (subject to additional constraints, in particular the CFL max 
constraint). 

In addition to the above criteria, based purely on a stability analysis, some other criteria 
for application of implicit schemes can be applied. For example, within boundary layers the 
implicit scheme may be preferred, because it provides faster convergence of the boundary 
fluxes and offers direct control of the natural boundary conditions. Many of these issues are 
yet to be studied. 

4.3 Artificial Dissipation 

In order to suppress spurious oscillations of the solution, a variety of models of artificial 
dissipation are used. In this work, we assume that the artificial dissipation can be introduced 
as the additional flux in the Navier-Stokes equations in the form: 

u + Ff^Fl + F*, 

where the artificial dissipation flux is the function of the solution vector and its derivatives: 

Ff = Ff{u,u tj ) 

with corresponding Jacobians defined as 

BF a 

pA gf « 

du 


The advantage of this approach is that the artificial dissipation can be treated using 
exactly the same formulation and procedures as for the natural viscosity. In the implicit 
algorithm, for the sake of generality, a fourth implicitness parameter 6 is introduced for the 
terms associated with the artificial dissipation. In the calculation of the stiffness matrices, 
right-hand sides and boundary terms, the same formulas are used as for the natural viscosity. 
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Within the above framework, various models of artificial dissipation can be formulated 
relatively easily. For example, a straightforward extension of the original Lapidus dissipation 
[26] to multidimensional cases yields artificial fluxes defined as 


Ft = k"u 


(4.23) 


with 

k = c k A e | 

where c k is a coefficient (usually between zero and one), A e is the element area, and v x are 
the components of velocity vector (no summation on i). The Jacobians P and R can be 
defined by a straightforward differentiation of formula (4.23). 

Another generalization of the Lapidus concept was proposed by Lohner, Morgan, and 
Peraire [27]. Within the framework proposed here, the fluxes corresponding to their model 
are of the form: Q 

or 

Ft = kliljUj ( 4 - 24 ) 

where l is the normalized gradient of the magnitude of velocity, 

j _ grad\v\ 

\grad\v\\ 

and the coefficient k is calculated from the formula 


k = c k A e {l ■ grad(v •/)) 


The Jacobians P and R can be defined by differentiation of formula (4.23). If, for simplicity, 
dependence of k and l on the solution is disregarded, then: 


Pi = 0 

Rij = kliljl 


(4.25) 


where / is the identity matrix of dimension M. In the incremental equation (4.11), the above 
approach leads to an additional term of the form 
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which differs slightly from the original formula proposed by Lohner and Morgan (equation 
(9) in [27]). These two versions are equivalent only for l constant throughout the domain. 
It can be verified, however, that the directional derivative used in [27] is not in divergence 
form and thus cannot be directly used in the variational formulation for arbitrary l. 

4.4 Implementation of Boundary Conditions 

Before discussion of boundary conditions for the implicit Taylor-Galerkin method, it is use- 
ful to quote the general result of Strikwerda [3], which specifies the number of boundary 
conditions necessary for well-posedness of the linearized Euler and Navier- Stokes equations. 
These results are summarized for two-dimensional problems in the table below. 


Type of 
Boundary 

Euler 

(not regularized) 

Navier-Stokes 

supersonic inflow 

4 ess 

4 ess 

subsonic inflow 

3 ess 

3 ess + 1 nat 

subsonic outflow 

1 ess 

1 ess + 2 nat 

supersonic outflow 

0 

0 ess + 3 nat 

no-flow 

1 ess 

1 ess 2 nat 

solid wall 



— isothermal 

— 

3 ess 

— heat flux 


2 ess + 1 nat 


In this table “ess” denotes the essential boundary conditions and “nat” denotes natural 
boundary conditions. The essential conditions are to be imposed on the characteristic vari- 
ables rather than the conservation variables. It is of importance to note that the numbers 
presented in the table are true for problems that are not regularized. If as is the cae of 
virtually all computational techniques — some artificial diffusion is built into the algorithm 
or added explicitly, natural boundary conditions should be imosed on these terms even for 
Euler problems. Moreover, since artificial diffusion (in contrast to the natural viscosity) 
affects all the conservation variables, the number of natural boundary conditions for these 
terms should actually be one more than for the (nonregularized) Navier-Stokes equations. 

In the jargon of finite difference methods, essential boundary conditions are equivalent to 
the boundary conditions to be specified, and the natural boundary conditions are equivalent 
to the boundary conditions to be extrapolated from the interior of the computational domain. 
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Within the finite element context, the basic idea of implementing these boundary conditions 
is based on the concept of constrained minimization. The essential boundary conditions 
are treated as constraint functions associated with the variational formulation. Penalty 
methods are usually applied to enforce these conditions in the original system and result in 
additional stiffness matrices and right hand load vectors. This method has been applied for 
implementing the open (inflow/outflow) boundary conditions, no-flow boundary conditions, 
and no-slip boundary conditions. 

In this section, the numerical implementation of pressure outflow boundary conditions 
and porous wall boundary conditions are presented. 

Enforcement of Prescribed Pressure 

For the case of subsonic outflow with a specified pressure there is one essential boundary 
condition to be specified, namely the prescribed pressure. The procedure currently being 
applied in this case is the following: 

(a) Apply a supersonic outflow procedure to account for corresponding boundary 
integrals (continuation from the interior condition). 

(b) Impose the one essential boundary condition (pressure) by the penalty method. 

The supersonic outflow procedure is described elsewhere and it amounts to rigorous calcu- 
lation of boundary integrals resulting from variational formulation. 

Beginning with the constitutive relation for the pressure 


p — (7 — 1 ) (E t — Ek) = 7 (E t — Ek) 


( 4 . 26 ) 


the condition for a prescribed pressure p is (in incremental form) 


A p = p- p n 


( 4 . 27 ) 
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The corresponding penalty term (in energy form) is 


2 l - [Ap - (p -p n )] 2 


(4.28) 


and in variational formulation is 


I [Ap - (p - p")] ■ S( Ap) 


(4.29) 


where the test function for pressure is its own variation. 

Now observe that Ap can be expressed in terms of conservation variables as 


dp 

Ap = 7 T- Au = d • Au 


(4.30) 


where d = vi,u 2 ,l} T (in the two-dimensional case). Therefore the penalty term 

(4.29) becomes 

l -[d- Au-(p-p n )](d-u) (4.31) 

where u is a test function. Using the standard requirement that the variational equation be 
satisfied for every v leads to the vector equation: 


- \d ■ Au — (p — p n )] d = 0 
£ 


(4.32) 
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and after regrouping 


-[<f (g> <£\ Au = 

'Cr> | 

a. 

1 

penalty stiffness matrix 

right hand side 


(4.33) 


Note: The first component of d has density in the denominator, so that this approach will 
blow up if the density is close to zero. To stabilize this procedure, one can multiply (4.33) 
by p 2 , to get 

(4.34) 



where d = pd = {£jt,mi,m 2 ,p}. 


Porous Wall Boundary Conditions 

Generally speaking, the burning surface of a solid propellant may be regarded as a porous wall 
with boundary motion. The burning rate and the mass flow rate across the burning surface 
are, theoretically, dependent on the local flow conditions such as pressure and temperature 
and the composition of the solid propellant. In order to implement the porous wall boundary 
condition with possible boundary motion simultaneously and consistently, a total number 
of four quantities are prescribed (see Fig. 4.5). These quantities are the tangential and 
normal velocities of the wall (urn' and u;vw), the mass flux across the wall (m/v), and the 
temperature of the wall (TV). The requirement that these prescribed quantities have to 
be satisfied at all times on the burning surface results in three boundary conditions to be 
imposed. 

(a) The tangential velocity of the fluid is equal to the tangential velocity of the wall, 

Uj = Ujw 


(b) The balance of the net mass flux, 

p(un — uhw) = fn/v = normal momentum 
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U NW 



Figure 4.5: Prescribed quantities on a burning boundary. 

(c) The temperature of the fluid is equal to the temperature of the wall, 

T = T W 


It is important to note that the normal velocity of the fluid is not equal to the normal 
velocity of the wall due to the prescribed mass flux. The difference of these two quantities 
contributes to the net mass flux across the wall. 

The numerical implementation of these boundary conditions results in the modifications 
of the stiffness matrices and global right hand side associated with boundary integrals. For 
condition (a), the incremental form consistent with the incremental time stepping algorithm 
can be expressed as 

Auj = uxw — Ut (4.35) 

The constraint function associated with this condition is defined by 


g = A ut — ( utw — u t) 


(4.36) 


The variational form of the penalty term resulting from the constraint function g (tested 
with its own variation Sg = SAuj) is 


- [Aur — ( utw — ^r)] ' &{&Ut) 
£ 


(4.37) 
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The function g tested with the variation of the incremental tangential velocity, (4.35), will 
result in penalty terms affecting both the continuity and momentum equations. To avoid 
this unphysical situation, the constraint function, g , may be tested with the variation of the 
tangential momentum, that is 


— [A uj — (utw — *4)3 ' ^Amy 


(4.38) 


where mj denotes the tangential component of momentum. Note that 


where 


A uj = 

dU T 

dU 

Aw = d/\u 

A my — 

dm j 
dU 

Aw = dAu 

duj 


ut T\ T2 

du 


i ■> •? 

p p p 

dm? 

du 

- [0, Ti, T 2 , 0] 


(4.39) 


(4.40) 


and (Ti,T 2 ) represent the x,y components of the unit tangential vector. Substituting (4.39) 
into (4.38), the penalty term becomes 


— [d ■ A u — ( utw ~ *4)3 ^ ' v 

£ 


(4.41) 


Since the global variational formulation of the problem (not shown here) has to be satisfied 
for any test function v , the corresponding kernel of the penalty term in the stiffness matrix 
can be found from (4.41) as 

(4.42) 


k = -d (g) d 

£ 


and the corresponding right hand side is 


r = 




(4.43) 
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The condition (b) to be imposed on the porous wall states that 

p (un — Utfw) = 


or 


mw — putfw = m N 


The incremental form of this condition can be written as 


(4.44) 


Amjv — A pufjw = m/v — ( m )v — p Uu Nw) 


(4.45) 


This results in a constraint function L, given by 


L = (Am at - A pu NW ) - [mN ~ ( m N ~ p n UNw)} = 0 (4.46) 


By applying a typical penalty procedure as described for imposing condition (a), we obtain 
( L is tested with 8m N ) the corresponding kernel, k, and the right hand side, r, as 


where 


-d ® d 

(4.47) 

£ 


- [m^/ — (m n N — p n UN\v)} d 

(4.48) 

— — = [— unw, Ni, N 2 , 0] 

oU 

d du = l °’ N " N 01 

(4.49) 


and (N u N 2 ) denote the x,y components of the outward unit normal vector. 

For condition (c), we apply the same penalty procedure to enforce the prescribed tern 
perature. The resulting penalty term is (tested with the variation of total energy) 


AT - (Tw -T n )] • 6(E t ) 


(4.50) 
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The temperature can be expressed in terms of total energy, E t , and kinetic energy, E k , as 


r = 7(7 - 0(7-7 


or 


t = 7(7-1) hr- . 2 


E t m x m x 

p 2? 


(4.51) 


dT 


where m, represents the momentum components. The standard vector d - g u can be 
derived from (4.51) as 

, dT _( E, EA 

d ' = a? = 7 r^ w 

dT _7 U , _ dT_ _ 

^ 2 dm-i p ’ 3 dm 2 p 

d 4 = ~ , 7 = 7(7 - 1) 

dE t p 

The corresponding kernel and right hand side result from this condition is 

k = - d ® d 

£ 


r = l -(T w -T n )d 


where 


d = 7 


E t E k 1_ 

P 2 P 2 P P P j 


d = [0, 0 , 0 , 1] 


5 Moving Grid/Eroding Boundary Algorithms 

As mentioned earlier, the surface of the solid propellant is treated as a porous wall charac- 
terized by mass injection and boundary motion. The mass injection rate and the speed of 
the boundary motion are dependent on the local flow conditions and the composition of the 
solid propellant. If combustion phenomena are not considered, these quantities are usually 
determined according to empirical data. 
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As combustion phenomena has not been included as part of this project, and empirical 
approach has been implemented whereby the motion /erosion of the solid propellant boundary 
is prescribed. The moving grid algorithm that incorporates the motions into the deformations 
of the domain can be stated as follows: 

For a given computational domain, ST, but <9f2i represent a part of the boundary which 
is subjected to a prescribed boundary motion and 80 2 represent the stationary component, 
such that 80 2 = dO i U 80 2 (see Fig. 5.1). 



Figure 5.1: Computational domain with a moving boundary. 

In order to prevent unacceptable mesh distortions near the moving boundary, we have 
implemented an algorithm whereby the whole computational mesh is stretched. The grid 
velocities u k for the stretching are calculated by solving a boundary-value problem defined 
by the Laplace operator (which has certain smoothing properties): 

VV = 0 in 0 , k = 1,2 

u n = u N on 80i 

uyv = 0 on 80 2 

where u k represents the k - th component of grid velocity and u N is a prescribed normal 

velocity. Note that only the normal component of velocity is prescribed on the moving 

and stationary boundaries. This allows the grid to stretch along the boundary and better 
adjust to the erosion process. However, corner nodes belonging to two different boundary 
sections have two linearly independent normal vectors and thus two components of velocity 
are automatically prescribed. Note that the two equations are coupled due to the specified 
normal boundary velocity, u 

A preconditioned block Jacobi-CG iterative method was employed to solve the final 
system of linear equations for grid velocity. Once the grid velocities for each node point 
are determined, the new nodal positions are updated by using the time step size and the 
calculated grid velocity. 

X k = X k + Dt x u k 
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Presently, two options for determining the time step size were implemented in: 

(a) In the first option, the time step size Dt is calculated from from the fluid stability 
criteria. Although an accurate transient solution may be obtained, the boundary 
motion may be extremely slow due to the small time step size and the relatively 
slow motion of the moving boundary. Thus, this version may require extremely 
long solution times and high computational costs to move the boundary large 
distances. 

(b) In the second option, we specify a certain amount of time to advance the moving 
boundary without solving the fluid problem. During the boundary motion, the 
grid velocities are recalculated periodically to account for changes of the normal 
vector n and to prevent mesh distortion. Currently, the criterion for the corre- 
sponding time interval Dt is that the mesh can move only a prescribed distance 
before recalculation of the grid velocities. This distance is usually taken as a 
fraction of the averaged element size. 


6 Adaptive Mesh Strategies and Data Management 
Schemes 

Adaptive methods are an efficient means for improving the accuracy of a computational 
solution. In the implementation of such methods, the relative accuracy of a solution is de- 
termined by calculating an error indicator for each element in some appropriate norm. Once 
an error is determined, adaptive methods are used to change the structure of the approxima- 
tion of the problem to reduce the error below a preassigned limit. Changing the structure of 
the approximation may involve increasing or decreasing the number of elements (h- methods), 
shifting the grid points (r-methods), altering the order of the local finite element approxi- 
mation (p-methods), or any of several other techniques. 

The h - adaptive scheme incorporated into the solid rocket booster code uses a normalized 
error indicator to estimate the relative magnitude of the local element error. This type of 
error indicator has proven to have the capability of capturing shock waves with variable 
strengths and shock wave interactions. This normalized error indicator is defined as 


0 e = h e )~ sup 
U »'=1,2 


dU 


dxi 


(6.1) 
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where h e is the measure of the element size and U is the independent variable such as density, 
pressure, etc., in the error calculations. 

Based on this local error indicator, the h-refinement and unrefinement strategy can be 
summarized by the following steps: 


An /i-Refinement/Unrefinement Method 


The /i-procedure involves the following steps: 


1. For a given domain fl, such as that shown in Fig. 6.1, a coarse finite element mesh 
is constructed which contains only a number of elements sufficient to model basic 
geometrical features of the flow domain. 

2. As the adaptive process is designed to handle groups of four elements at a time, a finer 
starting grid is generated by by a bisection process, indicated in Fig. 6.1b, to obtain 
an initial set of element groups. 

3. The numerical solution is calculated on this initial coarse grid, and the error indicators 
6 e are computed over all M elements in the grid. Let 


#MAX — 


max 6 e 

1 <e<M 


4. Next, the groups of elements are scanned and the group errors are computed 

p 

^group = S 

k=\ 

where e k is an element number in group k and P = 4. 

5. Error tolerances are defined by two real numbers, 0 < <*,/? < 1. If 


0 e > /3#max 

element e is refined. This is done by bisecting element e into four new subelements. 

If 

^GROUP — ^MAX 

group k is unrefined by replacing this group with a single new element with nodes 
coincident with the corner nodes of the group. 

This general process can be followed for any choice of an error indicator. Moreover, it 
can also be implemented with any prescribed frequency. 
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Cb) 


Figure 6.1: (a) Refinement and unrefinement of a four-element group, and (b) a coarse initial 
mesh consisting of a four-element group. 
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Data Structures 


An important consideration in all adaptive schemes is the data structure and associated 
algorithms needed to handle the changing number of elements, their node locations and 
numbers, and the element labels. 

As noted in the preceding paragraphs, the algorithm is designed to process (refine or 
unrefine) in groups of four elements at each local refinement/unrefinement step. Consider, 
for example, the case of an initial mesh of 20 square elements shown in Fig. 6.2. We assign 
to each element in this mesh an element number, NEL = 1,2,. . . ,NELEM and to each global 
node a lable NODE. This array, NODES(J,NEL) relates the local number J(J = 1,2, 3, 4) of 
element NEL to the global node number NODES. In addition, the coordinates Xj,Yj of 
each node are also provided relative to a fixed global coordinate system. These numbers are 
filed in two arrays: 

NODES(J,NEL) = the array of global node numbers assigned to node J of 

element NEL 

XCO(JCO,NODE) = the array of J CO — coordinates of global node NODE(JCO 

= 1 or 2) 

Suppose that an error indicator is computed that signals that an element should be 
refined, say element 11, in the example. We must have some system for assigning appropriate 
labels to the new elements and nodes. Toward this end, a convention has been established 
that defines the connectivity of the specified element with its neighbors in the mesh. This 
information is provided by a third connectivity array. 

Thus, the bookkeeping of element and node numbers evolving in a refinement process is 
monitored by the arrays NODES( . , . ) , XCO( . , . ) , NELCON( . , . and an array 
LEVEL(NEL) which assigns a level number to element NEL. Initially, the same level can 
be assigned to all elements, and this level is an arbitrary parameter prescribed in advance 
by the user. Thus, provisions are now in hand for an arbitrary dynamic renumbering of 
elements and nodes. If, for example, for the mesh in Fig. 6.2 (element 11) is to be refined, 
we proceed through the following steps: 

1. Loop over the neighbors of element 11 (which is made possible with the NELCON 
array), checking the level of the neighboring elements relative to the level of element 

li; 

2. If any neighboring element has a level lower than element 11, then element 11 cannot 
be refined at this stage; 

3. If element 11 can be refined (as is the case in Fig. 6.2), generate new element numbers 
(thus changing NELEM and new node numbers for unconstrained nodes); 
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Figure 6.2: Mesh, node, and connectivity numbering in a model problem. 
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4. Compute the connectivity matrix NELCON for the new elements; 

5. Adapt the connectivity matrices for the neighboring elements (since the refinement of 
element 11 has now changed this connectivity); and 

6. Interpolate the solution for the new nodes. 

Consider, as an example, the uniform grid of four elements shown in Fig. 6.3a and 
suppose that the error estimators dictate that element A is to be refined. Thus, A is divided 
into four elements, I, II, III, IV, as shown, and the solution values at the junction nodes, 
shown circled in the figure, are constrained to coincide with the averaged values between 
those marked X. Note that the connectivities change in the process, e.g., the connectivities 
4 and 8 of element B are different. 

Next, assume that an additional refinement is required, and that we must next refine 
element III. We impose the restriction that each element side can have no more than two 
elements connected to it. Thus, before III can be refined, element B must first be refined, as 
indicated in Fig. 6.3b. The constrained node B1 in Fig. 6.3a now becomes active, while node 
Cl remains a constrained node. With B bisected, we proceed to refine III into subelements 
a,/?, 7,(5, and new constrained nodes, again circled in Fig. 6.3c, are produced. In this case, 
only element B had to be refined first in order to refine III, but, in general, the number of 
elements that must be refined in order to refine a particular element cannot be specified. 


7 Turbulence Modeling 

Reynolds— Averaged Navier-Stokes Equations 

For most engineering analyses, only the mean motion of a turbulent flow is of interest. The 
governing equations of the mean motion of a turbulent flow are usually derived by applying 
the Reynolds decomposition technique and an averaging procedure to the Navier-Stokes 
equations. 

The time-dependent, mass-averaged, full Navier-Stokes equations can be expressed as 


dp 

dt 
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Figure 6.3: Sequence of refinements of uniform mesh. 
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where 


r,, = -p^-f(M + ^)^ + (f + ^)(|^ + S:) 

i. _ _ ( JL + ii- 'l — h is the specific enthalpy 

q ' ” \Pr PrJ dn' 

and u, is the turbulent eddy viscity, Pr and Pr, are the laminar and turbulent Prandtl 
numbers, respectively. These equations are exactly the same as for the laminar -case except 
as an additional eddy viscosity has been added to the molecular viscosity The eddy viscosity 
may be calculated using a wide range of turbulence models which vary from a ge raic mo e s 
to k-e models. We have selected a simple algebraic model [291 to implement within the 
context of the implicit/explicit flow solver described in Section 4. 

Algebraic turbulence models, although simple in formulation, are very difficult to imple- 
ment if complicated geometries are to be handled. The following sections present details 
of the numerical implementation of Prandtl-van Driest inner layer turbulence model in the 
context of adaptive unstructured grids. 

Prandtl-van Driest Turbulence Model 

In the Prandtl-van Driest turbulence model, the turbulent eddy viscosity is calculated ac- 
cording to: 

Hi = P* \w\*P 

where 

p = local fluid density 

|u)| = local magnitude of vorticity 

l — Prandtl mixing length 

and 

£ = kyD 

k = 0.4, von Karman constant 
y = distance normal to a solid (no-slip) wall 
D = van-Driest damping function 
= 1.0 - exp(-y + /A) 
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and 


A = 26.0, van-Driest constant 

y + = wall function = p w ■ v w y/p w 
p w = density on the wall 
p w = laminar viscosity on the wall 
v w = viscous wall velocity 
= y/fjw/ Pw) 

t w = shear stress on the wall 

It can be shown that if all flow variables are normalized by some reference quantities, for 
instance, 0oo for velocity, Poo for density, Lc for length, and for the viscosity, the formula 
for calculating fi w and v w will be modified as follows: 

m = Re oo p M P 

V w = Re oo ■ 7 " \ujpw 

For the case where a node has multiple length scales, say ML, the Prandtl mixing length 
will be calculated by taking the harmonic mean value of these length scales 

It should be mentioned that for structured grids used in finite difference methods, the data 
structures applied for indicating a wall are essentially built up by the grid indices with 
appropriate flags, for instance, / = 1, or, I = I MAX', J = 1, or J = JMAX. This simple 
data structure, together with an ad hoc assumption that the grid system is orthogonal, 
allows the calculation of the length scales associated with a grid point (or node) to be made 

relatively easy. 

In finite element methods, the calculations of p and \w\ is straightforward within an 
element. However, the calculation of Prandtl mixing length is very difficult, especially for 
adaptive unstructured grids, due to its strong dependence on the boundaries. This geometry 
depenence makes the calculation of length scales for a node very expensive. One reasonable 
method to reduce this dependence and to increase the computational efficiency is to build 
up a data structure that can be readily used to calculate the length scales for a given point. 
By having the data structure built up, the calculation of the Prandtl mixing length can be 
performed in a manner that is consistent with p and |u>| in an efficient manner. 

Requirements in Designing a Data Structure for Implementing an 
Inner Layer Turbulence Model 

In designing a data structure there are several requirements which should be kept in mind: 
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1. Efficient (ready to use) 

Once the data structure is built up, the length scales for a given point (nodal points 
or integration points) should be readily decoded from the data structure. 

2. Economic (minimum storage) 

The storage should be kept at a minimum and be flexible (dynamic allocation). 

3. Geometry independent (complex geometry) 

The data structure should be designed such that it is independent of the complexity 
of the geometry. 

4. Modularized (easy to couple with various flow algorithms) 

The data structure should be designed such that it could be easily coupled with any 
flow algorithms. 

5. Extendable, reusable (ready to be used with a two-equation model) 

The data structure should be designed such that it can be used with other turbulence 
models in which length scale is one of the key parameters for calculating turbulent 
eddy viscosity. 

6. General (multiple length scales) 

The data structure should be designed so that it allows for variable length scales for 
each node. 

7. Grid-structure dependent 

The data structure should be dependent on the global grid structure only. This implies 
that the data structure should be rebuilt only when the global grid structure has been 
changed, for instance, after a grid adaptation. 

8. Readable 

The arrays used in the data structure should be easy to read for engineers and the 
designer. 

9. Integrity, compactness, etc. 
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Criteria in Building Up the Data Structure 

It is important to set up geometrical criterion in building up the data structure for imple- 
menting the inner layer turbulence model. Without using these criteria, erroneous length 
scales and eddy viscosities may be determined. Presently, the following criteria are used. 

• A length scale for a node is a valid projection and/or the minimum distance to the 
viscous boundary (Fig. 7.1). 

• A valid length scale should not break the boundary of the computational domain (Fig. 
7.2). 

• The maximum allowable number of length scales per node is limited to four (4). 

• Each node has at least one length scale associated with a viscous boundary (Fig. 7.3). 

• If no valid projection can be found, the minimum distance to a wall will be taken as 
the length scale associated with a viscous boundary. 

• Length scales for each node are selected from all potential length scales that satisfy 
the above criteria starting from the minimum value. 

Design of a Data Structure for Implementing an Inner Layer 
Model 

With the above functional specifications and criteria in mind, we have designed a data 
structure for efficiently implementing an inner layer turbulence model. 

The first set of data is utilized to store information about the no-slip boundaries and the 
solutions on these boundaries. All the viscous boundaries are stored in a discrete form, that 
is, a boundary is composed by many line segments (this is true because bilinear elements are 
used in our finite element code). Each line segment will be referred to as a viscous panel. 
Each viscous panel can be identified by four integers, the element number it belongs to, the 
side number, and the node numbers (with respect to total number of viscous nodes) of its 
two end points. The reason for storing the node numbers of the end points for each viscous 
panel is that they could be readily used for calculating the wall function (j/ + ) for each interior 
node by using interpolation. 

The second set of data is utilized to store information that can be readily used to cal- 
culate the Prandtl mixing length for each node during the solution procedure. This set of 
data consists of a pointer, the number of length scales, and the values of the length scales 
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interior node 



We 7.1: A valid projection of an interior node with respect to a viscous boundary. 



Figure 7.3: The minimum distance is selected as the length scale. 
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associated with a node. For each length scale, the associated viscous panel number and the 
local coordinate of the projection are also stored. A pointer array was designed to take care 
of variable length scales, to efficiently allocate the storage and to direct access to the da a 
to be used to calculate the length scale. 

Details of the data structure are described as follows: 

COMMON /TURBZR/ RH0WALL(2*MXBCD) ,TAUWALL(2*MXBCD) , 

CMUWALL(2*MXBCD) , RTZDAT (2 ,2*MAXND) , 

TLENX(MAXND) 

COMMON /TURBZI/ NVISPAN.NVISNOD ,N0DVPAN(4,2*MXBCD) , 

NMVISP(MAXND) .ITZDAT (2 *MAXND) , 

ITZPTR(MAXND) 


DEFINITIONS OF ARRAYS 


New Flow Property: 

TLENX : TURBULENT LENGTH SCALE FOR EACH NODE 


Wall Variables: 

RHOWALL : DENSITY FOR EACH NODE ON VISCOUS BOUNDARY 

TAUWALL : WALL SHEAR STRESS FOR EACH NODE ON VISCOUS BOUNDARY 
CMUWALL : LAMINAR VISCOSITY FOR EACH NODE ON VISCOUS BOUNDARY 


Global Arrays: 

NMVISP : NUMBER OF VISCOUS PANELS ASSOCIATED WITH A NODE 

(NUMBER OF LENGTH SCALES FOR A NODE) , 

ITZPTR : THE POINTER FOR EACH NODE (IN ARRAYS RTZDAT AND ITZDAT) 

RTZDAT : VALUE OF Y-NORMAL AND LOCAL PARAMETRIC COORDINATES 

ON A VISCOUS PANEL FOR EACH NODE, 

ITZDAT : LIST OF VISCOUS PANELS ASSOCIATED WITH A NODE 
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Information About the Viscous Panel: 

NVISPAN : TOTAL NUMBER OF VISCOUS PANELS = NUMELD+NUMELE 
NVISNOD : TOTAL NUMBER OF VISCOUS NODES = NUMBCD+NUMBCE 
NODVPAN : INFORMATION TO BE DECODED TO IDENTIFY THE NODES ON 

EACH VISCOUS PANEL, 

1 : ELEMENT NUMBER, 2 : SIDE NUMBER 

3 : FIRST NODE NUMBER, 4 : SECOND NODE NUMBER 

Summary of Storage: 

1. A number of 7*MAXND words are required for storing information associated with 
the length scales for each node (assume averaged two length scales for each node). 

2. A number of 3*MXVISP words are required for storing solutions on the solid (no-slip) 
walls for efficiently interpolating solutions. 

Technical Comments 

For a grid system with MAXND=2500, it requires 140 kbytes of storage (assuming an average 
of two length scales for each node). This storage requirement is comparable to the storage 
required by the data structure designed by P. Rostand [30]: 90 kbytes for the same grid 
size and a “single” length scale for each node. However, it must be noted that Rostand’s 
data structure was designed for implementing both an inner layer and outer layer turbulence 
model and for simple geometries only (convex geometry was assumed). In addition, after a 
closer examination, one finds that there are some hidden storage requirements in Rostand’s 
calculations (for instance, solutions on the normals to be used for interpolation) not taken 
into account. This needed data requires additional storage and extra computational time if 
Rostand’s data structure were to be. generalized to handle more complicated geometries. We 
estimate that an additional 140 kbytes of storage may be required for our data structures if 
the outer layer turbulence model is to be implemented. 

Numerical Implementation 

The major tasks involved in the numerical implementation of the inner layer turbulence 
model are basically composed of two parts: 1) how to build up the data structure, and 
2) how to use the data structure to calculate the eddy viscosity. Flowchart 1 (Fig. 7.4) 
shows how the data structure is constructed in a step-by-step fashion. Flowchart 2 (Fig. 
7.5) presents a global view of how the inner layer turbulence model can be coupled with 
other flow algorithms with just a few extra subroutine calls. Once the data structure is built 
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from Pre-Processor 




Figure 7.4: Flowchart for constructing the data structures for inner layer turbulence model. 
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Flow Chart for Implementing Inner Layer Turbulence Model 



Figure 7.5: Flowchart for implementing inner layer turbulence model. 
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up, the calculation of the eddy viscosity for a given point (usually the integration point) is 
straightforward within the flow solver. 


8 Numerical Examples 

8.1 Supersonic Nozzle With Small Throat Radius of Curvature 

A 45°-15° convergent /divergent conical nozzle was used in this test case (see reference [31] 
for detailed nozzle specifications). This geometry is characterized by its rapid contraction in 
the convergent part, sharp wall curvature at the nozzle throat and near parallel at the nozzle 
exit. The experimental work for this type of nozzle flow has been performed by Cuffel, et al. 
[32] and numerical results have been reported by Serra [31]. 

An initial grid of 61 x 21 nodes was used to solve this problem. Flow was initialized by 
using the analytic solution concluded from a quasi-one-dimensional isentropic flow [33]. The 
nozzle was assumed to connect to a reservoir such that the flow condition at the nozzle inlet 
could be treated as a uniform subsonic inflow. This inlet flow condition remained constant 
during the solution process. At the nozzle exit, a supersonic outflow condition was specified. 
After 400 time steps (with a minimum cost option and CFLBOUND = 5.0), the flow pattern 
was well developed in the nozzle. The adaptive package was then switched on and the grid 
was adapted (every ten steps to the first level) for another 200 time steps with a = 0.50, 
0 = 0.70. For the next 100 time steps, the grid was adapted every 10 time steps to the 
second level with a = 0.40, /? = 0.60. Finally, the grid was adapted to the third level for 
another 100 time steps. The final adapted grid consists of 3930 elements and 3750 nodes. 

As evidenced by the adapted grid shown in Fig. 8.1a, the grid is automatically refined 
and is well aligned with the shock pattern in the divergent part of the nozzle. This pattern is 
best illustrated by the Mach number contours, as shown in Fig. 8.2c. An oblique shock was 
triggered by a junction between the circular circular-arc throat and the divergent section. 
This oblique shock hit the nozzle centerline and reflected back into the flow domain. The 
streamline distribution over the entire nozzle passage is shown in Fig. 8.1b. The rapid 
expansion of the flow in the region close to the nozzle throat can be seen from a closeup view 
of the Mach contours in Fig. 8.2d. The comparisons of experimental data [32], numerical 
results from Serra [31], and our prediction for the Mach number distribution along the nozzle 
centerline and the nozzle wall are presented in Fig. 8.3. 
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Figure 8.3: Comparison of Mach number distribution for a 45°-15° conical nozzle 



8.2 Viscous Flow Over a Sphere 

Flow over a sphere has been a popular benchmark problem for validating CFD codes due 
to the facts that the fluid physics for this problem is well understood [34] and that many 
numerical results are available for comparison [34,35]. The flow conditions selected for this 
case are M = 0.1 and Re = 100.0. The computational domain was first discretized by using 
a mixed structured/unstructured grid. After the flowfield developed to a certain stage, the 
grid was then adapted to the first level in the region close to the solid boundary and in 
the separation region. The surface of the sphere was treated as a no-slip isothermal wall. 
The axis of symmetry was treated as a no-flow or no-penetration boundary. Characteristic 
boundary conditions were imposed on the rest of the artificial boundaries. No artificial 
dissipation was added for this test case. 

Figure 8.4a shows the adapted grid which consists of 2457 elements and 2535 nodes. 
The recirculation region on the le^side of the sphere can be seen clearly from the velocity 
vector plot, as shown in Fig. 8.4d. From the streamline plot shown in Fig. 8.4c, it can 
be observed that the flow separates at an approximate angle of 123° (measured from the 
leading edge stagnation point) and that the dividing streamline extends into the wake region 
with the distance s/D = 0.8. The vorticity distribution and pressure distribution along 
the surface of the sphere (as shown in Fig. 8.5) agree extremely well with the data from 
[34,35]. It should be mentioned that, in the solid rocket booster code, the reference velocity 
used to nondimensionalize the Navier-Stokes equations is the speed of sound at farfield. The 
characteristic length is based on the diameter of the sphere. Therefore, the vorticity and 
the pressure coefficient calculated by the code must be rescaled according to the following 

formulae: 

1. for the vorticity 

U a 

id — 

2 Mqo 

where u> A is obtained from the code and u is the normalized vorticity using freestream 
velocity. 


1_ 

7 


where P A is obtained from the ADAPT2D™ code. 


2. for the pressure coefficient 


C„ = 


P Poo 

1 2 
ryp OO^OO 


M 2 

J oo 


Pa ~ 
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Figure 8 4- Flow over a sphere, M = 0.1, Re = 100.0. (a) Adapted structured/unstructured 
grid, (b) close-up view of adapted grid, (c) streamline plot shows separatum bubble, (d) 
close-up view of velocity vector plot in separation region. 
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Figure 8.5: Flow over a sphere, M = 0.1, Rt — 100. (a) vorticity distribution on the sphere 
surface, and (b) pressure distribution on the sphere surface. 


60 


8.3 Simulation of Vortex Shedding Due to Motor Inhibitor 

The purpose of this test case is to verify the axisymmetric capability of the code. The 
problem of vortex shedding due to a motor inhibitor protruding from the wall to the port 
flowfield studied by Majumdar, et al. [36] will be utilized as the benchmark problem. 

Due to the nature of axisymmetry, an axial/radial section of the SRM was taken as the 
computational domain. For the sake of fair comparison, the initial mesh size used for solving 
this problem is exactly the same as it was used in [36]. The flow conditions employed in the 
simulation are M = 0.09897, Re = 10 5 . No artificial dissipation was added for this case. 
Flow was initialized by using one-seventh power law velocity distribution. This type of 
velocity distribution was also used to specify the inflow condition. At the outflow, a pressure 
boundary condition was imposed. 

The grid was adapted to the first level during the solution process. The final grid consists 
of 722 elements and 772 nodes. As seen in Fig. 8.6a, grids were automatically refined along 
the no— slip boundary to resolve the viscous boundary layer. Strong circulation due to the 
inhibitor can be seen clearly from the velocity vector plot shown in Fig. 8.6b. The streamline 
plot indicates a large separation region formed in the lee-side of the inhibitor. The pressure 
distribution is displayed in Fig. 8.7a. These results qualitatively agree well with the results 
obtained by Majumdar, et al. [36] as exhibited in Fig. 8.8. 


8.4 Internal Flow in the Turnaround Duct of Space Shuttle Main 
Engine 

Due to the highly restricted space, the most important design specifications for the Space 
Shuttle Main Engine (SSME) are minimal weight and size. This requirement results in 
the complicated design of engine components which are usually combined with structural 
and geometrical complexities. A typical example is the 180° bend turnaround duct (TAD). 
The strong curvature of flow passage in the TAD will cause high levels of turbulence, flow 
unsteadiness and flow separation, etc. The ability to model turbulent flows is therefore im- 
portant for any CFD code in order to predict internal flow in the turnaround duct. However, 
as pointed out by Monson, et al. [29], turbulence closure models have been developed and 
optimized for external flows. The type of model needed and its importance for internal flows 
with strong curvature still remains to be established. The purpose of this test case is to 
verify the accuracy of the simple algebraic turbulence model that has been implemented 
within the code. 

Specification of the computational domain is given in [29]. Two test cases with Re = 10 
and Re = 10 6 have been performed. The initial grids for both cases are shown in Fig. 8.9. 
Flow was initialized by using the one-seventh power law velocity distribution [37]. The 
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Figure 8.7: Simulation of vortex shedding due to motor inhibitor, (a) pressure contours, and 
(b) streamlines. 
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inhibitor 



Figure 8.8: Simulation of vortex shedding to to motor inhibitor [36]. (a) velocity distribution 
and (b) pressure contours. 




maximum Mach number in this velocity profile is M = 0.1. This amounts to an averaged 
uniform flow at the inlet with M = 0.0875. This velocity is used to normalize the longitudinal 
velocity distribution. At the inlet of the TAD, the boundary condition was specified by using 
the one-seventh power law velocity distribution. The pressure outflow boundary condition 
was imposed at the exit of TAD. A fully implicit scheme was used to solve this problem. To 
let the flow develop smoothly, the CFL number was gradually increased from 1.0 to 50.0. 
After 400 time steps, the grid was adapted to the first level in the regions where active 
flow events are likely to occur, see Figs. 8.10, 8.11. The specified convergence tolerance 
was aligned within 600 time steps (the convergence tolerance is in the order of 10 6 ). The 
pressure contours, vorticity contours, velocity distribution and streamlines are shown in Figs. 
8.12, 8.13. The results shown here agree qualitatively very well with the numerical results 
predicted by Chen by using an extended k-e model as shown in Fig. 8.14. It is interesting to 
note that in the case of Re = 10 5 , a small secondary separation bubble was predicted by our 
code, which is consistent with the observation in the experiment performed by Sandborn (as 
shown in Fig. 8.15). The comparisons with the data predicted by other CFD codes [29] for 
the longitudinal velocity distributions at several axial locations are shown in Figs. 8.16-8.21. 
Good agreement can be observed for those locations far away from the separation bubble 
while some discrepancy can be seen in the regions close to the separation bubble. Further 
numerical studies are required to investigate this difference. 


8.5 Porous Cylinder with Nozzle (Planar Case) 

The internal flow in a nozzle is usually driven by the mass injection from the surface of 
the solid propellant and the positive pressure gradient between the rocket chamber and the 
environment. The purpose of this test case is to apply the porous wall boundary condi- 
tion and the pressure outflow boundary condition for simulating the internal flowfield in a 
cylindrical-port cold-flow model. This cylindrical-port model is connected to a nozzle which 
was constructed using piecewise analytic functions. This type of problem has been studied 
by Sabnis, et al. [38] and will be resolved here. 

The initial grid consists of 50 x 20 elements in the cylindrical part and 30 x 20 elements 
in the nozzle section. The flow conditions employed for this case are the mass injection 
rate = 0.0018 and Reynolds number = 8.0 x 10 5 . No artificial dissipation was added and 
laminar flow was assumed. Flow was initialized as a quiescent condition over the entire 
computational domain. The head end of the cylinder and nozzle wall were treated as a 
no-slip isothermal boundary. The axis of the cylinder was treated as a no-flow boundary. In 
order to drive the flow in the nozzle smoothly, the pressure at the outflow boundary was first 
set to a value slightly below its corresponding quasi-one-dimensional isentropic pressure. 
During the solution process, this value was gradually decreased until supersonic flow was 
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Figure 8.9: Initial grids used for simulating internal flow in the TAD. (a) Re - 10 5 , (b) 
Re = 10 6 . 


Figure 8.10: Adapted grid, Re = 10 s . 
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Figure 8.12: Internal flow in the TAD, Re = 10 5 . (a) velocity distribution, (b) streamlines, 
(c) pressure contours, and (d) vorticity contours. 
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Figure 8.13: Internal flow in the TAD, Re = 10 6 . a) velocity distribution, (b) streamlines, 
(c) pressure contours, and (d) vorticity contours. 
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Figure 8.14: Internal flow in the TAD (Y. S. Chen), Re — 10 s , k-e model, (a) velocity 
distribution, (b) streamlines, (c) pressure contours, and (d) vorticity contours. 
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Figure 8.15: Internal flow in the TAD (Sandborn), Re = 86,000.0, mean velocity profiles at 
different axial locations. 
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Figure 8.17: Normalized longitudinal velocity distribution at 6 — 90.0 deg 
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Figure 8.19: Normalized longitudinal velocity distribution at xj H - 2.0. 
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Figure 8.21: Normalized longitudinal velocity distribution at x/H = 12.0. 
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developed at the outflow. Once the flow becomes supersonic at the nozzle exit, the pressure 
outflow boundary condition was turned off. A method that is similar to the extrapolation 
boundary condition procedure used by the finite difference community is then imposed on the 
supersonic outflow boundary. Due to the nature of the subsonic flow within the cylindrical 
port field, the pressure disturbance will be propagated from the nozzle section upstream 
toward the head end and reflected back and forth within the chamber. After 1000 fully 
implicit time steps, this type of disturbance still can be observed in the numerical solution. 
The numerical solutions presented here were obtained after 2000 time steps. 

Figure 8.22 shows the grid and a plot of the streamlines. It can be seen clearly that 
the internal flow was driven by the mass injection from the surface of the cylinder. The 
comparisons of the normalized axial velocity distribution at several axial locations are shown 
in Figs. 8.23-8.26. Very good agreement between our code, MINT and experimental data 
are evident. 

The adaptive package was activated after 200 time steps and the grid was refined to the 
first level. The resulting mesh contained 2560 elements and 2594 nodes as shown in Fig. 
8.27. As can be seen from Fig. 8.27a, refinement is clustered in the near wall region in the 
divergent section of the nozzle to resolve the viscous boundary layer. The viscous boundary 
layer grows rapidly and significantly reduces the effective cross section area in the divergent 
part as shown in Fig. 8.27b. 


8.6 Porous Cylinder with Nozzle (Axisymmetric Case) 

To test the axisymmetric option in conjunction with a porous boundary condition, the 
previous test case was resolved with the grid size and all flow conditions, boundary conditions, 
and time stepping procedures exactly the same as those used in test case 8.5. 

Very similar numerical results were obtained as compared with the previous test case. 
Figure 8.28 shows the initial grid and the streamline. It can clearly be seen that the internal 
flow was driven by the mass injection from the surface of the cylinder. The pressure contours 
shown in Fig. 8.29 indicates that in the nozzle chamber the pressure is nearly constant and 
the pressure changes rapidly in the vicinity of the nozzle throat due to a fast flow expansion. 
Figure 8.30 shows the velocity distribution in the nozzle section. A relatively thin boundary 
layer extended from the nozzle wall can be observed. After the internal flow was developed to 
a certain stage, we turned on the adaptive options with the maximum level of grid refinement 
limited to the first level. The grid was adapted in the region close to the viscous boundary 
to resolve the rapid change of flow velocity as shown in Fig. 8.31. 
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Figure 8.22: Porous cylinder with a nozzle, (a) grid, (b) streamlines. 
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Figure 8.24: Normalized axial velocity distribution at Z/D = 6.64.. 
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Figure 8.27: Closeup view of (a) adapted grid, and (b) velocity profiles in the nozzle section. 
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Figure 8.28: Initial grid and the stream function contours. 
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Figure 8.29: Pressure contours. 
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Figure 8.30: Velocity distribution in the nozzle section. 
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8.7 Slotted Chamber with Nozzle (Axisymmetric Case) 

The purpose of this test case is to validate the axisymmetric modeling capability for simu- 
lating internal nozzle flows with mass injection. In addition, an algebraic turbulence model 
was applied. The nozzle geometry was constructed using piecewise analytic functions. The 
inhibitor was treated as a radial slot and attached to the cylinder (see [39] for detailed 
geometry specifications). 

The initial grid consists of 70 x 20 elements in the combustion chamber region and 30 x 20 
elements in the nozzle section. The flow conditions employed includes the normalized mass 
injection rate = 0.0027 and Reynolds number = 8.0 x 10 5 . The mass injection rate is 
held constant for all porous boundaries. The initial condition and the rest of the boundary 
conditions are exactly the same as the previous test case. The surfaces of the inhibitor are 
treated as porous boundaries. The numerical solutions presented here were obtained after 
4000 time steps. Figure 8.32 shows the grid an streamline plots. It can be seen that 

the internal flow was driven by the mass injection from the surface of the cylinder. A global 
velocity vector plot and details of the flowfield around the aft— end of the slot are displayed in 
Fig. 8.33a. The high speed flow injected from the slot to the chamber caused a reversed flow 
in the vicinity of the joint of the slot and the cylinder, see Fig. 8.33b. Chamber pressure 
is nearly constant and varies rapidly in the nozzle section as shown in Fig. 8.33c. The 
closed-up view of the velocity vector plots in the nozzle section and Mach number contours 
are shown in Fig. 8.34. These figures confirm that the large contraction of the convergent 
part leads to a rapid flow expansion. 

The comparisons of the normalized axial velocity distribution at several axial locations are 
shown in Figs. 8.35-8.43. Generally speaking, the velocity boundary layers at all locations 
predicted by our code are much thinner. This discrepancy may be caused by many reasons 
such as the grid resolution, distribution of the mass injection rate, the use of appropriate 
turbulence models, and flow unsteadiness, etc. Further numerical studies are necessary to 
explain these differences. 

After 4000 time steps, the grid was adapted to the first level. The final grid consists 
of 2315 elements and 2437 nodes. As expected, most of the adaptation occurred along the 
solid wall of the divergent part of the nozzle as shown in Fig. 8.44a. A very thin velocity 
boundary layer can be observed in Fig. 8.44b. 

8.8 Moving Grid Algorithm 

To test the moving grid algorithm, a 45°-15° convergent-divergent nozzle [31], with an 
extended cylindrical part for modeling the solid propellant section, was used. As mentioned 
in Section 5, if the grid motion is coupled with the fluids problem, the boundary motion could 
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Figure 8.35: Normalized axial velocity distribution at Z/D = 0.49. 
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Figure 8.36: Normalized axial velocity distribution at Z/D = 0.89. 
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Figure 8.38: Normalized axial velocity distribution at Z/D = 1.54. 
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Figure 8.39: Normalized axial velocity distribution at Z/Z) = 1.89. 



OF pOO'? 


yjw.f' 


98 







R/RH 


SHTTED CHAMBER KITH R MIZZLE. NBRMflLIZEO Mill VELIC1TY 
ORTA FAIR AFRFI-TR-06-KM 
NIRHALIZEO AXIRL VELICITV. Z/0=4.D8 


LEGEND 

O EXP. DfiTfl 
O HINT 
a ft0RPT20 


o 



Figure 8.41: Normalized axial velocity distribution at Zj D = 4.09. 
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Figure 8.43: Normalized axial velocity distribution at ZjD — 6.51. 
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Figure 8.44: Closeup view of (a) adapted grid, and (b) velocity profiles in the nozzle section. 
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be extremely slow and is not suitable for large boundary motions. Therefore, the moving 
grid algorithm, Option 2, was used. The initial grid is shown in Fig. 8.45. After 1.5 time 
units of advancement, the grid distribution is shown in Fig. 8.46. It should be noted that 
the Laplacian operator that governs the grid motion is a smoothing operator and, therefore, 
the grid motion does not affect the quality of the initial grid as can be seen from Fig. 8.46. 


9 Future Extensions 

The success of the two-dimensional and axisymmetric analysis package in modeling various 
benchmark problems as demonstrated in Section 8 suggests that extensions of current sim- 
ulation capability to realistic three-dimensional flows in cavities with eroding boundaries is 
feasible. In this section, we will briefly describe some possible extensions of this project with 
respect to the enhancement of the two-dimensional code and for the development of a fully 
three-dimensional capability for use in the design and analysis of solid rocket motors. 

Two-Dimensional Enhancements 

The two-dimensional code has numerous special capabilities for modeling subsonic to super- 
sonic flow regions. To enhance these current capabilities the following options are proposed. 

1. Develop a menu driven, x-window-based interactive preprocessor and postprocessor. 

2. Develop a user-friendly two-dimensional structured/unstructured mesh generation 
package. 

3. Optimize (vectorize/parallelize) the implicit/explicit solution module to provide peak 
performance on the MSFC Convex and/or Cray computers. 

4. Enhance the moving boundary algorithm to include a physically based regression ve- 
locity which is a function of the flow characteristics such as pressure and temperature. 

5. Implementation of higher— order turbulence models such as a k-i or k-£ model. 
Three-Dimensional Code 

Many of the configurations used in solid rocket motors are three-dimensional in nature due 
to the shape of the combustion chamber. Thus a three-dimensional analysis capability which 
extends the two-dimensional/axisymmetric capabilities would clearly be of great interest and 
benefit. The following steps outline the development of a three-dimensional Navier-Stokes 
solid rocket motor code: 
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Figure 8.45: Initial grid used for simulating mo 
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ing grid/eroding boundary algorithm. 




Figure 8.46: Grid distribution after 1.5 time units of advancement. 
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1. Develop three-dimensional discrete models of the Navier-Stokes equations for com- 
pressible flow in domains with boundaries unergoing arbitrary motions, particularly 
with burning or receding surfaces where the burning rate is determined by local flow 
physics. 

2. Develop a three-dimensional mesh generation package for moving boundary simula- 
tions. 

3. Develop a grid visualization package for viewing moving grids in three dimensions. 

4. Develop a complete three-dimensional anisotropic /i-adaptive package including the 
data structure, refinement /unrefinement package, and directional error estimates. 

5. Extend the current two-dimensional burning boundary algorithm to three dimensions. 

6. Assemble the functional packages developed in Tasks 1 to 5 into user-friendly three- 
dimensional analysis and design package for modeling solid rocket motors. 

7. Develop a set of code validation problems to be supported by experimental results 
supplied by MSFC (if available). 

8. Develop complete user documentation. 
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Appendices 

A Jacobians Due to Source Terms 


Four groups of additional Jacobian matrices are required in the axisymmetric formulation. 
1. Jacobian matrix due to inviscid source term, B , is defined as 


B = 


d&_ 
du 

= -B 

y 


B = 7 


0 0 0 0 

0 0 0 0 

-q 2 / 2 u v — 1 

0 0 0 0 

L. 

where 

7 = 7 — 1 , q 2 = u 2 + v 2 

and 7 is the ratio of specific heats. 

2. Jacobian matrix due to viscous source terms, T, is defined as 

t = ^ = -t = -(t 1 + -t 2 \ 

du y y \ y j 


where 


T x = 


0 0 


0 0 0 0 

A dpu 2 A pu dp 
p 1 dx p 3 dx 

A dpv 2\pv dp X dp A dp 

p 2 dy p 3 dy p 2 dx p 2 dy 


0 0 
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0 0 0 0 


t 2 


0 0 0 0 

2 Jtfl o ±i o 

p 2 p 


0 0 0 0 


3. Jacobian matrices due to viscous source terms, Q, are defined as 

~ dSV _ 1 o 

Qi — Tp: — — Q i 

oU'i y 

0 0 0 0 

0 0 0 0 


Qi = 


^ d 0 0 


0 0 0 0 
~ dSV l r> 

Q 2 — TjT — Q2 

on, 2 y 

0 0 0 0 


Q 2 


0000 

hn 0 — 0 

p 2 p 


0000 
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OF- 

4. The viscous Jacobians, P' = ~ 1 , are defined as 

du 

0 0 0 0 

-^21 ^22 "^23 ^ 


P 1 = 


dj\ 

du 


P 31 P32 P33 0 

P 4 \ P42 PI 3 P44 


where the matrix components of P 1 are given by 


pi 

■*21 


1 l ' x , o m l^.l , n^ m lpA 

— ~iiRm hl - Xm 2 , 2 + 2p/? 1- 2 A — — I 

P^ V P P ) 


pi 

•*22 


PH 

• — P,1 


p 1 

-'23 


A dp 
p 2 dy 


p 1 — _ 

-*31 — 


H . n P ,2 . 0 P ,1 , m \ 

— — m 2 1 — mi 2 + 2m! 1- 2m 2 1 

p 2 \ ' ' p p y j 


pi 

^32 


.P ( P *+1 

p v p y / 


pi 

- r 33 


P 

‘- JP.l 

P 2 


Pi 

"‘41 


uP 2 \ = uP 21 + uP 3 \ - (m,T n + m 2 r 21 ) 


+ 


P 2 C„ 


— (pe)i -(- 2uim lil -f- 2u 2 m 2i i + ^2e 3tXj 3 u 2 ^ p,i 


P 1 

-*42 


— + uP 22 + vP 32 + ( _m i.i + 2u iP,i) 

p P c v 


pi 
-* 43 


— —■ + uP 2 3 + UP 33 H 5 — (— m 2> ! + 2 u 2 p i i ) 

p P 


P 1 

^44 


■~o — P < 1 

P 2 Cv 


no 



Setting i = 2 in the definition above for the viscous Jacobians gives: 


p> = ?h 

du 


0 0 0 0 
P 2 \ PI Ph 0 

Pi, Pi PI o 

p2 p2 p2 p2 

r 4\ r 42 r 43 r 44 


where the matrix components are defined by 
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B Methods for Treating Constrained or Hanging 
Nodes 

Suppose that one is interested in performing element calculations for the element NEL shown 
shaded in Fig. B.l. The solution within this element can be expressed as 


U e = £ tiNODES(/.NEL)ty/ 


;= i 


(B.l) 


= Wj#! + U 2 tf 2 + «3®3 + «4®4 


where $/ is the local shape function associated with node I , and u { is the solution vector at 
node 1. However, notice that nodes 1 and 3 are hanging nodes and the numerical solution 
at these two nodes are constrained by 


Ul = 


Uz 


U 4 + Ug 

2 

U 4 + U8 


I 


(B.2) 


Substituting (B.2) into (B.l), we have 


u 


= + "?) 4l + u 2 $2 + ( B -3) 


Rearranging terms, we obtain 


u 


= Ug -y + tt 2 ^2 + “By + “ 4 (y + Y + ^ 4 ) 


(B.4) 
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Figure B.l: Elements consisting of constrained nodes. 


By defining new shape functions, 'F', as 




/ 

1 


Hi 

2 


*2 = *2 



(B.5) 


and defining (u 9 , u 2 , « 8 , u 4 ) as the physical degrees of freedom associated with this element, 
then the solution vector within NEL can be interpreted as 


4 

U e = ^UNODES (I,NEL)^I 

1=1 


(B.6) 


where 

N0DES(1,NEL) = 9 , N0DES(2,NEL) = 2 

NODES(3,NEL) = 8 , N0DES(4,NEL) = 4 

and Vj are the modified shape function defined by (B.5). The method being applied to 
handle the constraint nodes for the Cartesian problems is based on equations (B.5) and 
(B.6). That is, one modifies the data structure such that element NEL was defined by the 
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physical degrees of freedom (nodes 9, 2, 8, 4) rather than its original geometrical degrees of 
freedom (nodes 1, 2, 3, 4). The local shape functions for these physical degrees of freedom 
have to be consistently modified according to equation (B.5). 

In the axisymmetric case, the conservation variable u = yu was interpolated separately 
for y and u. It is important to note that y is a purely geometrical quantity. The use of the 
modified shape functions to interpolate y within an element will cause a severe inconsistency 
for the conservation variable u. Therefore, y should be interpolated by using the original 
element shape function, that is 

y e = Yly 

: = 1 

Thus, equations (B.5), (B.6), and (B.7) completely define the method for treating the con- 
strained nodes for the axisymmetric problem. 


C Projection of Surface Nodes 

The basic idea behind an adaptive methodology is that during the solution process, the grid 
is automatically refined and/or unrefined according to certain error criteria. These criteria 
are developed such that when the solution process is repeated using an updated grid the 
accuracy of the solution is improved. 

Due to the refinement and unrefinement of the grid, many new nodes and elements 
are generated and some existing nodes are eliminated. This in general does not cause any 
difficulties if the new nodes are generated interior to the domain but may lead to problems if 
the nodes are adjacent to a prescribed boundary and no precautions are taken to ensure that 
the new grid points lie on the given surface. Obviously, it is necessary that as the refinement 
continues the resolution of the profile of the object must be improved as well as the accuracy 
of the solution. Without a precise interpolation or projection, these new nodes generated on 
the prescribed boundaries during the adaptive procedures can destroy the accuracy of the 
original geometrical representeation of a two-dimensional object and thus the solution. 

If the surface of a two-dimensional object can be expressed by a few simple algebraic 
equations, an accurate interpolation or projection is easy to perform. Unfortunately, for 
most realistic problems the surfaces of a two-dimensional object are usually defined by a 
finite number of discrete points rather than a set of continuous algebraic functions. Thus 
we are faced with the problem, how does one accurately and efficiently project nodal points 
onto the surface of the two-dimensional object defined by a finite number of discrete points. 
One popular method used by the CAD/CAM industry [40] which we will also employ is 
described below. 
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The Parametric Cubic (PC) Curve in Space 

The algebraic form of a PC space curve can be expressed as 

3 

r(u ) - 5 >u‘ 
i=o 


where r = (x,y,z) T , u is the parameter that characterizes the curve and we assume 
from 0 to 1 , and a, are the constant coefficients. 

Suppose that the end conditions t are given, that is 

r( 0 ) = a 0 

r(l) = a 0 + ai + 02 + 03 ^ 

r'(0) = oj 

r'(l) = 3a 3 + 2a 2 + a! 


or in matrix notation 


1 

0 

1 


O 

O 

O 

h— * 

l 



r(l) 


1111 


a 2 

r'( 0 ) 


0 0 10 


a i 

, r'O) . 


3 2 1 0 _ 


do 


Solving (C.3) for a, yields 


^3 


' r ( 0 ) 

a 2 

a a 

= (M) 

r(l) 

r'( 0 ) 

do 


. r'(l) _ 


where 


M = 


2-211 
-3 3 -2 -1 

0 0 10 
10 0 0 


Substituting (C. 4 ) into (C.l) and rearranging, we obtain the geometric form of r 


r(u) = r(0)6i(ti) + r(l)fra(«) + r'( O)fts(ti) + r'(l) 6 4 («) 


(C.l) 

varies 


(C.2) 


(C.3) 


(C.4) 


(C.5) 


(C. 6 ) 
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where 


bi(u) = 2 u 3 - 3u 2 + 1 ' 
b 2 (u) = —2u 3 + 3u 2 > 

63(11) = u 3 — 2u 2 + u 
6 4 (u) = u 3 -u 2 


(C.7) 


are the so-called blending functions. In matrix notation, equation (C.6) can be written as 


r(u) 



(C.6) 


Replacing r by x. y, and z in (C.8), the full representation of a PC space curve becomes 


x(u) 

y( u ) 
z{u) J 


u 3 u 2 ul 


(M) 


x(0) 2/(0) 2(0) 

x(l) y{l) z(l) 
x'(0) y'(0) x'(O) 

Lx'(l) y # (l) ^(1) J 


(C.9) 


Note that the PC plane curve is the special case of (C.9) by neglecting the ^-coordinate. 


Segmenting a PC Curve 

In most practical engineering applications, PC curves are usually defined in a piecewise 
manner. Assuming that the end values of a segmented PC curve are given by u x and u 2 , 
then the geometric coefficients of the segmented curve in terms of the given curve are modified 


■ H(0) ' 


r(uj) 

K(l) 


r(u 2 ) 

«(0) 


Aur'(ui) 

fi'(l) 


Aur'(it 2 ) 


where Au = u 2 - m and R denotes the 1, y, and 2 coordinates for the segmented PC curve, 


Projection of Surface Nodes 

Usually, the arc length coordinate is commonly used as the parameter that characterizes a 
PC curve. This implies that those points to be supplied to construct a PC curve have been 
pre-sorted in certain order such that every point on the curve can be assigned a unique arc 
length coordinate. The element connectivity array as explained in Section 6 can be used to 
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identify if a node lies on the prescribed surface. For a point to be projected onto a PC curve, 
the method we applied is first to find the closest segment in the PC curve to the point. Once 
this segment was identified, the projection of the surface node is equivalent to the following 

statement: 

For a given point with Cartesian coordinate r = (x,y), find its arc length coordinate, s, 
and associated Cartesian coordinate r(s) such that 


d = ||r(s) - r| 


(C.ll) 


is a minimum, where 


r(s) 


s 2 si] M 


As = s2 


r(sl) 


r(s2) 

(C.12) 

Asr(sl) 

Asr(s2) _ 


— si 

(C.13) 


and M is given by equation (C.5). 

Obviously, to minimize d is equivalent to find the zero-root of a fifth order polynomial 
f(s) defined as 

f(s) = [r{s)-r]-r'{s) (C.14) 

Newton- Ralphson’s iterative method may be applied to determine its zero-root as follows: 


1. Use s = 0.5(sl + s2) as an initial guess. 

2. Calculate f{s) and /'(s), where f'(s) = \r(s) - r] • r"(s) + r'(s ) • r’(s) 

3. Calculate 8 = — — \ 

f{s) 

4. Update s = s + 8 

5. Check if |<5| < tolerance, 
if yes, root = s and stop 
if no, go to step 2. 

By knowing the arc length coordinate of an arbitrary point on the PC curve, its corre 
sponding Cartesian coordinate can be easily determined by using equation (C.12). 


117 



D Interface With GAMMA2D 

Presently, the interface with the SRB2D code is through a user supplied grid file which con- 
tains node definition, element definition, element connectivity and, possibly, some geometric 
data for defining user prescribed profiles such as a nozzle contour or the profile of an eroding 
pocket. These profiles are treated as parametric cubic splines so that ruled curves may be 
readily established to perform surface nodal projections (for grid refinement) and surface 
node redistribution (for remeshing), see Appendix C for details. COMCO has developed 
an in-house grid generator, GAMMA2D, to automatically generate this grid file which can 
be readily interfaced with SRB2D. For details see the GAMMA2D user’s manual and the 

SRB2D user’s manual. 


E Summary of the Postprocessing Capabilities 

For the sake of completeness of the final report, we summarize the current postprocessing 
capabilities of the code. 

1. Pointwise or nodewise aerodynamic data extraction. 

For the purpose of engineering analysis, it is highly desirable that the CFD codes can 
provide quantitative information about the flowfield. Nodewise data extraction allows 
a user to extract desired flow prooperties from the computed solution for a given point 
or a given boundary in the computational domain. These aerodynamic data include 
all conservation variables, all primitive variables, pressure coefficient, Mach number, 
shear stresses, entropy change, total enthalpy change, total energy loss, etc. 

2. Image of the initial grid and adapted grid with an option to show a section of the 
computational domain.. 

3. Contour plots. 

Contour plotting of flow variables is a very popular method to show the distribution of 
flow properties throughout the computational domain. The code is capable of plotting 
the contours for the following flow properties: all conservation variables, all primitive 
variables, Mach number, vorticity, entropy change, total enthalpy change, total energy 
loss, molecular viscosity, eddy viscosity, turbulent length scale, etc. 

4. Velocity vector plot and streamline plot. 

These plots are very useful to resolve complicated flow topologies such as vortex shed- 
ding, flow separation, flow reattachment, etc. 
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5. Total mass flux, total force, total heat transfer, total moment calculations. 

To activate these postprocessing capabilities in the code, the user provides key words in 
semantic form in the input deck. These key words are usually accompanied with some user 
specified options. See the SRB2D user’s manual for details. 
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